Role of bond orientational order in the crystallization of hard spheres 
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With computer simulations of the hard sphere model, we examine in detail the microscopic path- 
way connecting the metastable melt to the emergence of crystalline clusters. In particular we will 
show that the nucleation of the solid phase does not follow a two-step mechanism, where crystals 
form inside dense precursor regions. On the contrary, we will show that nucleation is driven by fluc- 
tuations of orientational order, and not by the density fluctuations. By considering the development 
of the pair-excess entropy inside crystalline nuclei, we confirm that orientational order precedes po- 
sitional order. These results are at odd with the idea of a two-step nucleation mechanism for fluids 
without a metastable liquid-liquid phase separation. Our study suggests the pivotal role of bond 
orientational ordering in triggering crystal nucleation. 



I. INTRODUCTION 

After almost 150 years from the pioneering works of 
Gibbs, a satisfactory understanding of the crystalliza- 
tion process is still lacking. Classical Nucleation The- 
ory (CNT) describes nucleation as an activated process 
in which the free energy gain of the solid phase over the 
fluid one is in competition with the free energy penalty of 
forming an interface for nuclei of size smaller than a char- 
acteristic length, the critical nucleus size. The two main 
limitations of CNT have been recognized as the following: 
i) the capillarity approximation, i.e. the assumption that 
the thermodynamic properties of the small crystal nuclei 
are the same as the bulk crystal, such as specific volume 
or surface tension; ii) all the order parameters involved 
in the transition proceed simultaneously, and the transi- 
tion can be described effectively by one order parameter 
(as density for the nucleation of liquid droplets from a 
supersaturated gas). Many efforts have been done in or- 
der to relax some of the assumptions of CNT, showing in 
fact that nuclei do not form at bulk condition and that 
the process is non-classical, i.e. the order parameters in- 
volved in the transition do not change simultaneously [l[ . 

On these grounds, the two-step nucleation mechanism 
has emerged as a candidate description for many of the 
discrepancies between CNT and experiments. The idea is 
the following. The transition from the fluid phase to the 
solid phase involves at least two order parameters, one 
associated with the breaking of translational symmetry 
and the other with the breaking of rotational symme- 
try. In the crystallization of spherical particles, usually 
these order parameters are identified with density and 
bond orientational order. In the two-step mechanism sce- 
nario, the first fluctuation to trigger the transition is den- 
sity fluctuation, which leads to the formation of a dense 
droplet in the melt. Then the droplet starts to be struc- 
tured and a crystal is formed. Large-amplitude density 
fluctuations can emerge from the criticality of a nearby 
liquid-liquid (L-L) phase transition, which is the case of 
globular proteins [2( or colloids with small range attrac- 
tions Q. This also led to the speculation that critical 
fluctuations could enhance the crystallization rate [3, Q . 
Rather surprisingly, the same mechanism was claimed 
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Recently we have provided a detailed microscopic 
study of the nucleation pathway in ultrasoft potential 
systems [l(| and in hard sphere systems [ll|, both hav- 
ing purely repulsive potentials and thus without a L- 
L transition. These works highlighted the role played 
by bond orientational order in the crystallization pro- 
cess |12|, finding no evidence of the dense precursors 



upon which the two-step assumption is built. On the con- 
trary, the transition resembles a two-step process where 
the first step is the formation of extended structured 
regions of high orientational order, which progressively 
densify. The bulk density is in fact reached only at a 
very late stage, when the nucleus becomes several times 
the critical nucleus size. The first step, structuring at 
(almost) constant density is achieved due to the wet- 
ting of bond-oricntational regions to small crystalline nu- 
clei. The identification of the bond orientational order 
as the appropriate coordinate to describe the nucleation 
stage has unveiled that the polymorph selection stage 
starts already in the metastable fluid phase before nu- 
cleation occurs [l(| [H| • With bond orientational order 
we have also studied the competition between crystalline 
structures and icosahedral (or five-fold symmetric) struc- 
tures [Ullill, i.e. the mechanism by which crystallization 
is prevented and that paves the way to out-of-equilibrium 
states (i.e. glasses). 



In the present contribution we will focus on the nu- 
cleation pathway in monodispcrsc hard spheres. We 
will show that, while density fluctuations have a very 
short lengthscale, bond orientational order has an ex- 
tended lengthscale that increases with supersaturation. 
The nucleation events are localized inside the regions of 
high bond orientational order. The process is thus non- 
classical, with orientational order preceding translational 
order as the nuclei grow. 
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II. METHODS 

We run Umbrella Sampling Monte Carlo simula- 
tions of N = 4000 monodispcrse hard spheres in the 
isothcrmal-isobaric (NpT) ensemble. Hard spheres are 
ideal for studying crystallization and have already pro- 
vided tremendous contributions to our basic understand- 
ing of crystal nucleation [l4T - [l8j . In the following, 
lengths are given in unit of the diameter a, and pressure 
in units of k^T/a 3 , where fcs is the Boltzmann constant. 
We place the spheres randomly in a simulation box at 
packing fraction n = 0.5352 and equilibrate the system 
at reduced pressure dpa 3 = 17.0. At this pressure the 
liquid is metastable with respect to crystallization, with 
a difference in chemical potential between the liquid and 
solid state of 8\A(i\ = 0.54 pj. We run Umbrella Sam- 
pling simulations, where a biasing potential is added to 
the system Hamiltonian to sample crystalline clusters of 
different sizes, in order to extract information on the free 
energy barrier and the critical cluster size. Details on the 
implementation can be found in Ref. [20j . Configurations 
for clusters at different sizes are extracted from the simu- 
lations and analyzed. The Umbrella Sampling biasing is 
used to gain better statistics by stabilizing nuclei at the 
desired size. But we confirm that the same results are 
obtained if the configurations are extracted from direct 
simulations (still possible at this pressure) without bias- 
ing potential. At 0pa 3 — 17 the free energy barrier is 
0AF ~ 18 and the size of the critical nucleus is n c ~ 80. 

The identification of crystalline particles follows the 
usual procedure [2l|. A particle is identified as crystal if 
its orientational order is coherent (in symmetry and in 
orientation) with that of its neighbors. The £-fold sym- 
metry of a neighborhood around each particle i is charac- 
terized by a (21 +1) dimensional complex vector (q;) as 



J2f=i' ^m(fy), where I is a free integer 



parameter, and m is an integer that runs from m = —£ 
to m = £. The functions Yi m are the spherical harmonics 
and f y is the vector from particle i to particle j. The 
sum goes over all neighboring particles Nb(i) of parti- 
cle i. Usually Nb(i) is defined by all particles within a 
cutoff distance, but in an inhomogeneous system the cut- 
off distance would have to change according to the local 
density. Instead we fix Nb(i) = 12 which is the num- 
ber of nearest neighbors in icosahedra and close packed 
crystals (like HCP and FCC) which are known to be the 
only relevant crystalline structures for hard spheres. The 
scalar product (q&(i)/\q&(i)\) ■ (q&(j)/\q&(j)\) quantifies 
the similarity of the two environments. If it exceeds 0.7 
between two neighbors, they are deemed connected. We 
then identify a particle as crystalline if it is connected 
with at least 7 neighbors plj . 

The liquid-to-solid transition is characterized by the 
symmetry breaking of both translational and orienta- 
tional order 22[. Translational order expresses the rel- 
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in the top row of Fig. [T] A good order parameter for 
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FIG. 1: Translational order (upper row) and orientational 
order (bottom row) in a fictitious liquid (left column) and 
crystal (right column) configurations. Translational order ex- 
presses the relative spacing between particles in the system, 
while orientational order the relative orientation between the 
neighbors around each particle. 



translational order is the local density p, which is calcu- 
lated by constructing the Voronoi diagram and measur- 
ing the volume of the cell associated with each particle. 
More generally, translational order can be obtained from 
two-body correlation functions, such as the pair correla- 
tion function g(r) |22l [. Since we are dealing with hard 
sphere systems, where there is no energy involved, the 
most meaningful definition of translational order comes 
from the two-body excess entropy, defined as: 



§2 



-| / dr[g(r)log(g(r)) - g(r) + 1} 



(1) 



Orientational order, instead, expresses the relative ori- 
entation between the neighbors around each particle, as 
shown in the bottom row of Fig. [T] Unlike transla- 
tional order, which is obtained from two-body correla- 
tion functions, orientational order is obtained by consid- 
ering many-body correlations. For hard spheres it is well 
known that the local bond-order analysis introduced by 
Steinhardt [23j is an adequate measure of orientational 
order. This order parameter is obtained by construct- 
ing a rotational invariant from the quantity qe m (i) pre- 
viously introduced. In the present work we are going to 
focus on the coarse-grained quantity Qe, which can be 
obtained in the following way. First qe m (i) is spatially 
coarse-grained [2~H 
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where N(i) are the neighbors of particle i. Then the 
following invariant is obtained 



1 



47T 

2Z + 1 



(3) 



3 




0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 



<?6 



le-01 1 10 

FIG. 2: Probability distribution for the metastable phase in 
the p — Qe space. The dashed black line is a contour line. The 
dashed white arrow is instead a steepest descent path from the 
maximum to a high Qe point of the probability distribution 
function. 

where we focus on the case I = 6 (corresponding to the 
six-fold symmetry of the crystal structures). The effect of 
coarse graining Q 6 has proven effective not only in reduc- 
ing the signal-to-noise ratio [Hj] > but also eliminates the 
signal from five- fold symmetric structures [HI |2o| which 
do not add coherently. In this way the order parameter 
goes continuously from low values in the fluid phase to 
high values in the crystal phase. 

III. RESULTS 

We start by looking at the two-dimensional probabil- 
ity distribution of density p and bond oricntational order 
Qe, for a metastable fluid state at pressure f3pa 3 = 17 
(before the appearance of the critical nucleus), reported 
in Fig. [2] The probability distribution is related to the 
Landau free energy, F(Qe,p) = — 1:bT log P(Qe, p)- The 
free energy can be well fitted with a full cubic polyno- 
mial, for which the most important term is of the form 
Qep 2 - This term is responsible for the shape of contours 
lines (dashed line in Fig. because the interaction is 
linear in Qe and quadratic in p, the system can increase 
its orientational order without an increase of its transla- 
tional order, but the contrary is not true, and an increase 
in density also increases the average Qe- A similar term 
was proposed in order to describe quasi-crystals forma- 
tion, where the ordering of the orientational field does 
not imply an ordering of the translational field [2(| . Note 
also that a small linear coupling between Qe and p exists 
at high Qe, seen for example in the small slope of the 
steepest descent path (white dashed arrow) in Fig. [2] As 
we will discuss later, we believe this small coupling term 
to be responsible for previous observations of two-step 
nucleation processes in hard spheres. 

Next we address the correlation length of these order 
parameters. In the two-step scenario, crystal nucleation 
occurs inside dense droplets that form because of fluc- 
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FIG. 3: Same configuration but with different sets of particles 
plotted, (left) 5% of particles having the highest value of p. 
Colors are given by sorting all the particles (also the ones not 
plotted) according to their value of Qe, and dividing the set in 
two halves: red (dark gray) particles have a high value of Qe, 
while green (light gray) particles a low value of Qe- (right) 
5% of particles having the highest value of Qe- Colors are 
given according to the value of p of each particle: red (dark 
gray) particles have a high value of p, while green (light gray) 
particles have a low value of p. 



tuations in the density field. In a recent study we have 
compared correlation lengths for both translational and 
orientational order [l3| and found that, while the corre- 
lation length associated with density (and any two-body 
quantity) is always very short and does not grow with 
increasing supcrsaturation, the correlation length of ori- 
cntational order is an increasing function of supersatura- 
tion. This already suggests that the mechanism of forma- 
tion of dense droplets is unlikely in hard spheres, while 
a viable mechanism would be the nucleation inside local- 
ized regions of high orientational order [ll|; E3 ■ We show 
this in Fig. [3J by simply looking at a typical configuration 
in the metastable state. On the left panel we plot a sub- 
set of particles (5% of the total) having the highest value 
of p. It is immediately clear that the position of parti- 
cles is approximately uncorrelated, without any apparent 
lcngthscale. The colors are given according to the bond 
orientational order field, with high Qe particles being col- 
ored in red (dark gray), and low Qe particles colored in 
green (light gray) . We see that there is a lack of correla- 
tion of particles in dense environments and the value of 
Qe- Going back to Fig. [5] this is confirmed by observing 
that, at high density, the Qe value (while increasing on 
average) is equally likely to have low or high values. On 
the right panel of Fig.[3]we plot instead the subset of par- 
ticles (5% of the total) having the highest value of Qe- 
Unlike the density field, bond orientational order clearly 
displays a strong correlation between the particles, with 
structures resembling droplets. In this case, the colors 
are given according to the value of p of each particle: red 
(dark gray) particles having a high value of p and green 
(light gray) particles a low value of p. This time there 
is a clear correlation between high values of Qe and high 
values of p, so clusters of high oricntational order will 
appear denser. This correlation comes from the linear 
term that we noted in the distribution function of Fig. [T] 
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FIG. 4: Radial profile of the averaged critical nucleus (nuclei 
with size between 80 and 90 particles). The profiles are for 
different order parameter (OP) fields (from top to bottom): 
OP = Qa (black line), OP = p (red or dark gray line) and 
OP — S2 (blue or light gray line). The fields are normalized 
as to be in the fluid phase and 1 in the bulk crystal phase. 
As usual, distances are in unit of the hard sphere diameter a. 



As we will show soon, since crystals are born from fluc- 
tuations of the orientational order field, a measure of the 
density of the regions from which the crystal originates 
would give a value higher than the average melt density. 
But this does not mean that the crystal originated due 
to a fluctuation of the density field, but it is simply a 
consequence of this small linear coupling between p and 
Qe- 

We now proceed to study the nucleation stage, where 
a crystal is nucleated up to the critical size, eventually 
overcoming the free energy barrier and starting the crys- 
tallization process. We first focus on nuclei of a size com- 
parable to the critical nucleus size (so we average over 
many independent configurations where the biggest nu- 
cleus is of size between 80 and 90 particles). In Fig. [4] 
we plot the radial profile for three different fields, Qq, 
p and the two-body excess entropy S2, as a function of 
the distance from the center of the nucleus (/). The 
fields are normalized as to be in the fluid phase and 
1 in the bulk crystal phase. First we note that the pro- 
cess is non-classical, where all fields have still not reached 
their bulk value, and the radial distribution is different 
for each field. In particular, going from the fluid phase 
(high value of r') to the center of the crystal nucleus 
(r J = 0), the Qq field is the first to develop, while density 
is lagging behind the development of the orientational 
field. This again proves that nucleation is driven by ori- 
entational order and not by density. The same result is 
confirmed by looking at other definitions of translational 
order, which can be extracted from two-body correlation 
functions. For hard spheres systems, where there is no 
energy term in the free energy, a physically motivated 
definition of translational order is given by the two-body 
excess entropy S2 defined in Eq. [TJ To obtain the radial 
dependence of S2(r') we introduce the functions g(r;r') 
which measure the radial distribution function averaged 
over particles at distance r' from the critical nucleus. In 



o Mi tzl -I 

1 1.5 2 2.5 3 

r 

FIG. 5: Pair correlation function g(r;r') calculated at dif- 
ferent distances r' (in steps of Ar' = 0.5<r) from the critical 
nucleus center. The g(r, r') at the center of the critical nucleus 
(r = 0) is the one closest to the dashed line, which is the pair 
correlation function for a bulk rhep crystal. Distances are in 
unit of the hard sphere diameter a. 



formula 

- N N 

pg{r;r') = -( Yl X^( r - r i + r i)>* 

t' <ri<r' +dr' j=^i 

where the ensemble average (■•■)* is over Umbrella Sam- 
pling configurations with a crystal nucleus of critical size. 
The two-body excess entropy function is then 

S2(r') = -| J dr[g(r;r')log(g(r;r'))-g(r;r') + l] 

where the integral is carried up to a distance correspond- 
ing to the second nearest-neighbors shell (as nuclei are 
still small). The result shown in Fig. 2] clearly indicates 
that the translational order (expressed by S2) is lagging 
behind orientational order (expressed by Qe). 

The slow positional ordering inside the critical nucleus 
can be seen in Fig. O where the pair correlation function 
g(r, r') is plotted at different distances (r') from the cen- 
ter of the critical nucleus. The g(r,r') at the center is 
the one closest to the dashed line, which is the pair cor- 
relation function g{r) of a bulk rhep crystal. Even at the 
center of the critical nucleus (r' = 0), the pair correlation 
function is far from that of a bulk crystal. 

Already from Fig. [4] it is clear that the two-step sce- 
nario with dense precursors is unlikely, as orientational 
order precedes the densification process. To have a di- 
rect confirmation of the absence of a two-step process, in 
Fig. [6] we plot the density at the center of the nucleus as 
a function of the nucleus size. The increase in density is 
linear with the size of the nucleus, and extrapolates to 
the fluid density in the limit of vanishing size of the nu- 
cleus. This result is in contrast to the expectation of the 
two-step nucleation mechanism, which instead predicts 
a non-linear dependency of the density on the nucleus 
size: first an enrichment at constant size, followed by a 
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FIG. 6: Excess density at the center of the nucleus (Ap = 
p(r = 0) — p fluid) as a function of the size of the nucleus 
R, with R c — 2.3 a being the critical nucleus radius. The 
increase in density from the fluid phase up to nuclei bigger 
than the critical size is linear, with Ap(R = 0) = 0. 



size growth with little density change. Figurc|6]thus con- 
firms that there is nothing anomalous with the density 
increase during the nucleation process. 



IV. CONCLUSIONS 

In the present article we have examined in detail the 
possibility of a two-step nucleation process from the su- 
percooled state in hard spheres. The two-step nucleation 
is based on a specific behavior of two order parameters: 
density and structure. The fluid phase first forms a dense 
precursor, which is then structured to form a crystalline 
nucleus. Nucleation is then initiated by density fluctua- 
tions, and followed by the superposition of a structural 
fluctuation. Unlike this scenario, we stress that hard 
spheres take another route to crystal nucleation. Nu- 
cleation is promoted by extensive fluctuations in bond 
orientational order. These fluctuations act as precursors 
for the formation of crystalline nuclei. The increase of 
density inside these regions lags behind the development 
of orientational order, so that one can effectively have 
small crystalline nuclei with densities still very far from 
bulk conditions. The reason why the transition from a 
fluid structure to a solid structure occurs with little den- 



sity change is due to the fact that nuclei form in regions 
of high orientational order, and so are wetted by fluid 
regions with a high value of this field. The structuring of 
the nuclei then requires just small adjustments of particle 
positions, without big density changes fllj ] . 

This decoupling between the orientational order pa- 
rameter and the translational order parameter has also 
important consequences for the glass transition. Once 
the positional order transition is avoided, orientational 
order keeps growing with increasing supersaturation. It 
has been recently shown that this increase of orienta- 
tional order is linked with the slowing down of the dy- 
namics, displaying critical-like behavior (25l. \2l\. 

Interestingly, the mechanism by which the positional 
ordering transition is avoided seems to be also correlated 
with the development of orientational order. We have re- 
cently shown that crystalline structures are not the 
only structures that are stabilized by an increase of ori- 
entational order. There is in fact a competition between 
crystalline structures and icosahedral structures (five- 
fold symmetric structures), and we have shown that both 
an increase in density and polydispersity [l3| stabilize the 
icosahedral structure, effectively suppressing crystalliza- 
tion. This again confirms the effect of frustration on crys- 
tallization [28| , as in the case of 2D spin liquids 123 . |30| , 
and as recently observed in metallic glasses [3lL |32| | - 

We stress that the two-step scenario is instead fully 
consistent for systems which display a L-L transition. In 
this case, the presence of a minima (even if metastable) in 
the system's free energy corresponding to the liquid phase 
can start the nucleati on p rocess without a macroscopic 
phase separation [2hfl l33l |34|. 
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